/* _NSPROJ_DPP_C_
 * 
 * SOLVES THREE DIMENSIONAL, INCOMPRESSIBLE NAVIER-STOKES
 * EQUATION WITH RECTANGULAR DOMAIN USING FINITE DIFFERENCING 
 * ON A STAGGERED GRID
 * UTILIZES THE PROJECTION METHOD.  THE CONSERVATIVE FORM
 * CONVECTION IS DISCRETIZED WITHOUT THE USE OF UPWIND TREATMENT; 
 * BOTH THE ADVECTION AND DIFUSSION TERMS ARE TREATED EXPLICITLY. 
 * THE CURRENT SETUP IS FOR A LID DRIVEN CAVITY PROBLEM WITH A WELL
 * KNOWN VELOCITY AT THE TOP AND NO-SLIP BOUNDARY CONDIITIONS AT
 * ALL OTHER BOUNDARIES
 * 
 * ORIGINALLY WRITTEN BY DR. ZHI-GANG FENG: UNIVERSITY OF NORTH TEXAS
 * 
 * PORTED TO C AND 3D BY JOHN LEIDEL: GRADUATE STUDENT
 * 
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <math.h>
#include "nsproj.h"
#include "mpi.h"


//-----------------------------------------START_SIM
// START_SIM
//--------------------------------------------------
int start_sim( char *datafile, int argc, char **argv )
{
    /* VARS */
    double Re          = 0;     //-- Reynolds Number
    double dt          = 0;     //-- Time Step
    double ct          = 0;     //-- Current Time
    double tf          = 0;     //-- Final Time
    double lx          = 0;     //-- Width of Box
    double ly          = 0;     //-- Height of Box
    double lz          = 0;     //-- Depth of Box
    double dx          = 0;     //-- delta_x
    double dy          = 0;     //-- delta_y
    double dz          = 0;     //-- delta_z
    double simtime     = 0;     //-- TOTAL SIMULATION WALLCLOCK TIME
    double lsimtime    = 0;     //-- TOTAL ITERATION WALLCLOCK TIME
    double step1       = 0;     //-- TOTAL TIME FOR STEP_1
    double step2       = 0;     //-- TOTAL TIME FOR STEP_2
    double step3       = 0;     //-- TOTAL TIME FOR STEP_3
    double iot         = 0;     //-- TOTAL TIME SPENT IN I/O
    double bt          = 0;     //-- TOTAL TIME SPENT IN BOUNDARY CONDITIONS
    double cmt         = 0;     //-- TOTAL TIME SPENT IN COMMUNICATION
    int nt             = 0;     //-- NUMBER OF ITERATIONS [TOTAL_TIME/dT]   
    int c              = 0;     //-- Counter
    int size           = 0;     //-- MPI_COMM size
    int rank           = 0;     //-- MPI_COMM rank
    int rtn            = 0;     //-- MPI_RTN
    int tec            = 0;     //-- DETERMINES TEC OUTPUT
    int nsp            = 0;     //-- DETERMINES NSP OUTPUT
    int prov           = 0;     //-- PROVIDED

    struct timeval starttime;   //-- TIMER
    struct timeval endtime;     //-- TIMER
    NSMESH mesh;                //-- MESH [ALL CELLS]
    NSMESH lmesh;               //-- LOCAL_MESH [ONLY LOCAL CELLS]
    char project[1024];         //-- PROJECT NAME 

    NSPROCS procs;              //-- PROCESSOR BREAKOUT
    /* ---- */

    /* INIT MPI ENV */
    rtn = MPI_Init_thread( &argc, &argv, MPI_THREAD_MULTIPLE, &prov );
    if( rtn != MPI_SUCCESS )
    {
	fprintf( stderr, "%s\n", "ERROR : COULD  NOT INITIATE MPI ENVIRONMENT" );
	return -1;
    }

    MPI_Comm_size( MPI_COMM_WORLD, &size );
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    /* SANITY CHECK THE PROCESSOR COUNT */
    if( size > _MAX_MPI_THREADS_ )
    {
	if( rank == 0 )
	{
	    printf( "%s%d\n", "ERROR : EXCEEDED MAXIMUM NUMBER OF SCALABLE PROCESSORS: ", size );
	}
	MPI_Finalize();
	return -1;
    }

    if( rank == 0 )
    {
	printf( "%s\n", "NSPROJ_DPP v. 0.2" );
    }
    
    /* PARSE CONFIG FILE */
    if( rank == 0 )
    {
	if( nsproj_readdatafile( datafile, 
				 &Re, 
				 &dt, 
				 &tf, 
				 &lx, 
				 &ly, 
				 &lz, 
				 &mesh.uN,
				 &mesh.uS,
				 &mesh.vE,
				 &mesh.vW,
				 &mesh.wF,
				 &mesh.wB,
				 &mesh.mx, 
				 &mesh.my, 
				 &mesh.mz,
				 &nsp, 
				 &tec,
				 project ) != 0 )
	{
	    printf( "%s%s\n", "ERROR : COULD NOT READ DATA FILE : ", datafile );
	    
	    return -1;
	}
	
	/* SETUP GRID POINTS && TIMING */
	nt = (int)ceil( tf/dt);
	dt = tf/nt;
	dx = lx/mesh.mx;
	dy = ly/mesh.my;
	dz = lz/mesh.mz;

	mesh.nx = mesh.mx+2;
	mesh.ny = mesh.my+2;
	mesh.nz = mesh.mz+2;
	
	/* SANITY CHECK VARIABLES */
	if( (dt<=0) || (dx<=0) || (dy<=0) || (dz<=0) || (lz<=0) || (ly<=0) || (lz<=0) || (nt<1) )
	{
	    printf( "%s\n", "ERROR : DATA OUT OF BOUNDS" );
	    
	    return -1;
	}
    }

    MPI_Barrier( MPI_COMM_WORLD );

    if( rank == 0 )
    {
	printf( "%s\n", "...DETERMINING MESH DISTRIBUTION" );
    }

    /* DETERMINE GRID DISTRIBUTION */
    if( nsmpi_getprocdistro( &procs, mesh, size, rank ) != 0 )
    {
	if( rank == 0 )
	{
	    printf( "%s\n", "ERROR : COULD NOT DETERMINE PROCESSOR/MESH DISTRIBUTION" );
	}
	MPI_Finalize();
	return -1;
    }

    MPI_Barrier( MPI_COMM_WORLD );

    if( rank == 0 )
    {
	printf( "%s\n", "...DISTRIBUTING RUN DATA TO ALL RANKS" );
    }

    /* DISTRIBUTE INITIAL DATA */
    if( nsmpi_distributedata( &mesh, 
			      &Re, 
			      &dt, 
			      &tf, 
			      &nt,
			      &lx, 
			      &ly, 
			      &lz,
			      &dx, 
			      &dy, 
			      &dz, 
			      rank ) != 0 )
    {
	if( rank == 0 )
	{
	    printf( "%s\n", "ERROR : COULD NOT DISTRIBUTE DATA TO OTHER THREADS ");
	}
	MPI_Finalize();
	return -1;
    }

    MPI_Barrier( MPI_COMM_WORLD );

    if( rank == 0 )
    {
	printf( "%s\n", "...ALLOCATING FULL MESH" );
    }

    if( rank == 0 )
    {
	/* ALLOCATE FULL MESH */
	if( nsproj_meshalloc( &mesh ) != 0 )
	{
	    if( rank == 0 )
	    {
		printf( "%s\n", "ERROR : COULD NOT ALLOCATE MEMORY FOR MESH" );
	    }
	    
	    MPI_Finalize();
	    
	    return -1;
	}
    
	
	/* INITIATE MESH */
	if( nsproj_meshinit( &mesh ) != 0 )
	{
	    if( rank == 0 )
	    {
		printf( "%s\n", "ERROR : COULD NOT INITIALIZE MESH ELEMENTS" );
		nsproj_meshfree( &mesh );
	    }
	    
	    MPI_Finalize();
	    
	    return -1;
	}
    }
    
    MPI_Barrier( MPI_COMM_WORLD );

    if( rank == 0 )
    {
	printf( "%s\n", "...ALLOCATING LOCAL MESHES" );
    }
    
    /* ALLOCATE LOCAL MESHES */
    if( nsmpi_alloclocalmesh( &lmesh, procs, rank ) != 0 )
    {
	if( rank == 0 )
	{
	    printf( "%s\n", "ERROR : COULD NOT ALLOCATE MEMORY FOR LOCAL MESH ELEMENTS" );
	    nsproj_meshfree( &mesh );
	}

	MPI_Finalize();
	return -1;
    }
    
    /* BEGIN SIMULATION */
    if( rank == 0 )
    {
	printf( "%s\n", "....BEGINNING SIMULATION" );
	printf( "%s\n", "....TIMING IN SECONDS" );
	printf( "\n%s\n", "ITER\t CUR_T\t\t STEP_1_T\t STEP_2_T\t STEP_3_T\t BC_T\t\t COM_T\t\tIO_T\t\t ITER_T\t\t TOTAL_T" );
    }
    
    MPI_Barrier( MPI_COMM_WORLD );
    
    for( c=0; c<nt; c++ )
    {   
	/* DISTRIBUTE DATA TO ALL RANKS */
	gettimeofday( &starttime, NULL );
	nsmpi_distributemesh( &mesh, &lmesh, procs, rank, size );
	MPI_Barrier( MPI_COMM_WORLD );
	gettimeofday( &endtime, NULL );

	cmt += (double)(endtime.tv_sec-starttime.tv_sec);
        cmt += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	/* SOLVING TENTATIVE USTAR/VSTAR/WSTAR  WITH IMPLICIT METHOD */
        gettimeofday( &starttime, NULL );
        velstar( Re, dt, dx, dy, dz, &lmesh );
	MPI_Barrier( MPI_COMM_WORLD );
       	gettimeofday( &endtime, NULL );
	
	step1 += (double)(endtime.tv_sec-starttime.tv_sec);
        step1 += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	/* SOLVE FOR PRESSURE USING SOR */
	gettimeofday( &starttime, NULL );
	pressure( dt, dx, dy, dz, &lmesh );
	MPI_Barrier( MPI_COMM_WORLD );
	gettimeofday( &endtime, NULL );

	step2 += (double)(endtime.tv_sec-starttime.tv_sec);
        step2 += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	/* PRESSURE CORRECITON */
	gettimeofday( &starttime, NULL );
	pressure_correction( dt, dx, dy, dz, &lmesh );
	MPI_Barrier( MPI_COMM_WORLD );
	gettimeofday( &endtime, NULL );

	step3 += (double)(endtime.tv_sec-starttime.tv_sec);
        step3 += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;
		
	/* COLLECT ALL VALUES */
	gettimeofday( &starttime, NULL );
	nsmpi_collectmesh( &mesh, &lmesh, procs, rank, size );
	MPI_Barrier( MPI_COMM_WORLD );
	gettimeofday( &endtime, NULL );
	
	cmt += (double)(endtime.tv_sec-starttime.tv_sec);
        cmt += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	/* BOUNDARY CONDITIONS */
	gettimeofday( &starttime, NULL );
	if( rank == 0 )
	{
	    boundary_cond( dt, dx, dy, dz, &mesh );
	}
	MPI_Barrier( MPI_COMM_WORLD );
	gettimeofday( &endtime, NULL );
	
	bt += (double)(endtime.tv_sec-starttime.tv_sec);
        bt += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;
	
	/* OUTPUT ALL CURRENT DATA */
	gettimeofday( &starttime, NULL );
	
	if( rank == 0 )
	{
	    /* CURRENTLY ONLY WRITES DATA FROM RANK=0 */
	    if( nsp )
	    {
		nsproj_writedata( project, mesh, ct );
	    }

	    if( tec )
	    {
		nsproj_writetecdata( project, mesh, c );
	    }
	}
	gettimeofday( &endtime, NULL );

	MPI_Barrier( MPI_COMM_WORLD );

	iot += (double)(endtime.tv_sec-starttime.tv_sec);
	iot += (double)(endtime.tv_usec-starttime.tv_usec)*0.000001;

	lsimtime = step1 + step2 + step3 + iot + bt + cmt;
	simtime += lsimtime;

	/* PRINT STATUS TO SCREEN */
	if( rank == 0 )
	{
	    printf( "%d%s%f%s%f%s%f%s%f%s%f%s%f%s%f%s%f%s%f\n", 
		    c, "\t",ct, "\t",step1, "\t",step2, "\t",step3, "\t",bt, "\t",cmt, "\t", iot, "\t", lsimtime, "\t",simtime );
	}

	/* RESET LOCAL TIMING */
	step1    = 0;
	step2    = 0;
	step3    = 0;
	lsimtime = 0;
	iot      = 0;
	bt       = 0;
	cmt      = 0;
	ct       = ct+dt; /* UPDATE CURRENT TIME [t] */
	
	MPI_Barrier( MPI_COMM_WORLD );
    }

    if( rank == 0 )
    {
	printf( "%s%s%s\n", "NSPROJ :: ", project, " COMPLETE!" );
	printf( "%s%f\n", "TOTAL SIMULATION WALLCLOCK TIME = ", simtime );
    }

    nsproj_meshfree( &lmesh );

    /* FREE MESHES */
    if( rank == 0 )
    {
	nsproj_meshfree( &mesh );
    }
    
    MPI_Barrier( MPI_COMM_WORLD );

    MPI_Finalize();

    return 0;
}

//-----------------------------------------MAIN
// MAIN
//---------------------------------------------
int main( int argc, char **argv )
{
    /* VARS */
    char *datafile = NULL;
    int rtn        = 0;
    /* ---- */

    if( argc < 2 )
    {
	fprintf( stderr, "%s\n", "USAGE: <mpirun> nsproj_mpp <datafile>" );
	return -1;
    }

    datafile = (char *)malloc( 1024 );
    if( datafile == NULL )
    {
	fprintf( stderr, "%s\n", "ERROR :  COULD NOT ALLOCATE MEMORY FOR DATA FILE" );
	return -1;
    }

    sprintf( datafile, "%s", argv[1] );
    
    rtn = start_sim( datafile, argc, argv );
    
    free( datafile );
    datafile = NULL;

    return rtn;
}
